home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_12_06 / gage / divdif.c < prev   
Encoding:
C/C++ Source or Header  |  1994-02-06  |  5.0 KB  |  191 lines
  1. /* Divided Difference Interpolation */
  2. /* Copyright 1994 by Philip Gage */
  3.  
  4. #include <stdio.h>
  5. #include <math.h>
  6.  
  7. #define POINTS 10     /* Max number of data points */
  8.  
  9. /* Function prototypes */
  10. void interpolate(int,int,int,double*,double*,double*);
  11. double evaluate (int, int, double*, double);
  12. void print_function (int, int, double*);
  13. void print_table(int,int,int,double*,double*,double*);
  14.  
  15. /* Generate polynomial or exponential function
  16.      polynomial    1=polynomial, 0=exponential
  17.      degree        Degree of desired solution
  18.      n             Number of data points
  19.      x,y           Data point values
  20.      coefficient   Output coefficient array */
  21.  
  22. void interpolate (int polynomial, int degree, int n,
  23.            double *x, double *y, double *coefficients)
  24. {
  25.   double xpower, data[POINTS], diff[POINTS];
  26.   int i, j, k;
  27.  
  28.   /* Copy input y values to data array */
  29.   for (i = 0; i < n; i++)
  30.     data[i] = y[i];
  31.  
  32.   /* Initialize coefficients to identity element */
  33.   for (i = 0; i <= degree; i++)
  34.     if (polynomial)
  35.       coefficients[i] = 0.0;
  36.     else
  37.       coefficients[i] = 1.0;
  38.  
  39.   /* Solve each high order term */
  40.   for (i = degree; i >= 0; i--) {
  41.  
  42.     /* Copy data values to difference array */
  43.     for (j = 0; j < n; j++)
  44.       diff[j] = data[j];
  45.  
  46.     /* Compute divided differences or quotients */
  47.     for (j = 1; j <= i; j++)
  48.       for (k = 0; k < n-j; k++)
  49.         if (polynomial)
  50.           diff[k] = (diff[k+1] - diff[k]) / 
  51.                     (x[k+j] - x[k]);
  52.         else
  53.           diff[k] = pow(diff[k+1] / diff[k],
  54.                         1.0/(x[k+j] - x[k]));
  55.  
  56.     /* Average current difference row */
  57.     coefficients[i] = 0.0; 
  58.     for (j = 0; j < n-i; j++)
  59.       coefficients[i] += diff[j]; 
  60.     coefficients[i] /= n-i; 
  61.  
  62.     /* Remove high order term from data */
  63.     for (j = 0; j < n; j++) {
  64.       /* Compute x[j] raised to i power */
  65.       xpower = 1.0;
  66.       for (k = 0; k < i; k++)
  67.         xpower *= x[j];
  68.       if (polynomial)
  69.         data[j] -= coefficients[i] * xpower;
  70.       else
  71.         data[j] /= pow(coefficients[i],xpower);
  72.     }
  73.   }
  74. }
  75.  
  76. /* Evaluate function at xvalue using Horner's rule */
  77. double evaluate (int polynomial, int degree, 
  78.                  double *coefficients, double xvalue)
  79. {
  80.   double yvalue;
  81.   int i;
  82.  
  83.   if (polynomial)
  84.     yvalue = 0.0;
  85.   else
  86.     yvalue = 1.0;
  87.   for (i = degree; i >= 0; i--)
  88.     if (polynomial)
  89.       yvalue = yvalue * xvalue + coefficients[i];
  90.     else
  91.       yvalue = pow(yvalue,xvalue) * coefficients[i];
  92.   return yvalue;
  93. }
  94.  
  95. /* Print symbolic function terms using ^ for power */
  96. void print_function (int polynomial, int degree, 
  97.                      double *coefficients)
  98. {
  99.   int i, identity, printing = 0;
  100.   char operator1, operator2;
  101.  
  102.   /* Setup for polynomial or exponential */
  103.   if (polynomial) {
  104.     identity = 0.0;
  105.     operator1 = '+';
  106.     operator2 = '*';
  107.   }
  108.   else {
  109.     identity = 1.0;
  110.     operator1 = '*';
  111.     operator2 = '^';
  112.   }
  113.  
  114.   printf("\nF(X) = ");
  115.   for (i = degree; i >= 0; i--) {
  116.     if (coefficients[i] != identity) {
  117.       if (printing)
  118.         printf(" %c ",operator1);
  119.       printf("%g",coefficients[i]);
  120.       if (i > 0)
  121.         printf("%cX",operator2);
  122.       if (i > 1)
  123.         printf("^%d",i);
  124.       printing = 1;
  125.     }
  126.   }
  127.   printf("\n");
  128. }
  129.  
  130. /* Print table of results */
  131. void print_table (int polynomial, int degree, int n,
  132.           double *x, double *y, double *coefficients)
  133. {
  134.   int i;
  135.   double result, error, percent;
  136.  
  137.   printf("%15s%15s%15s%15s%15s %\n",
  138.          "X","Y","F(X)","Error","Error");
  139.   for (i = 0; i < n; i++) {
  140.     result = evaluate(polynomial,degree,
  141.                       coefficients,x[i]);
  142.     error = result - y[i];
  143.     if (y[i] != 0.0)
  144.       percent = 100.0 * error / y[i];
  145.     printf("%15g%15g%15g%15g%15g %\n\n",
  146.            x[i],y[i],result,error,percent);
  147.   }
  148. }
  149.  
  150. void main()
  151. {
  152.   int i, n, degree, polynomial, choice;
  153.   double x[POINTS], y[POINTS], coefficients[POINTS];
  154.  
  155.   /* Read number of data points from 2 to POINTS */
  156.   do {
  157.     printf("Number of data points: ");
  158.     scanf("%d",&n);
  159.   } while (n<2 || n>POINTS);
  160.  
  161.   /* Read each data point from user */
  162.   for (i = 0; i < n; i++) {
  163.     printf("x%d: ",i);
  164.     scanf("%lf",&x[i]);
  165.     printf("y%d: ",i);
  166.     scanf("%lf",&y[i]);
  167.   }
  168.  
  169.   /* Generate solutions until Quit selected */
  170.   while (1) {
  171.     printf("1 Polynomial interpolation\n");
  172.     printf("2 Exponential interpolation\n");
  173.     printf("3 Quit program\n");
  174.     printf("Choice: ");
  175.     scanf("%d",&choice);
  176.     if (choice < 1 || choice > 2) break;
  177.     polynomial = (choice == 1);
  178.  
  179.     /* Read degree of solution from 1 to n-1 */
  180.     do {
  181.       printf("Degree: ");
  182.       scanf("%d",°ree);
  183.     } while (degree < 0 || degree >= n);
  184.  
  185.     /* Compute function and print results */
  186.     interpolate(polynomial,degree,n,x,y,coefficients);
  187.     print_function(polynomial, degree, coefficients);
  188.     print_table(polynomial,degree,n,x,y,coefficients);
  189.   }
  190. }
  191.